﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "qrrt")]
    public static unsafe int qrrt(IntPtr a_ptr, int n, IntPtr xr_ptr, IntPtr xi_ptr, double eps)
    {
        double* a = (double*)a_ptr.ToPointer();
        double* xr = (double*)xr_ptr.ToPointer();
        double* xi = (double*)xi_ptr.ToPointer();

        return qrrt(a, n, xr, xi, eps);
    }

    /// <summary>
    /// 多项式方程求根QR方法
    /// </summary>
    /// <param name="a">a[n+1]存放n次多项式的n+1个系数。</param>
    /// <param name="n">多项式次数。</param>
    /// <param name="xr">xr[n]返回n个根的实部。</param>
    /// <param name="xi">xi[n]返回n个根的虚部。</param>
    /// <param name="eps">eps控制精度要求。</param>
    /// <returns>函数返回在求上H矩阵特征值时返回的标志值。若>0则正常。</returns>
    public static unsafe int qrrt(double* a, int n, double* xr, double* xi, double eps)
    {
        int i, j;
        double* q = stackalloc double[n * n];

        for (j = 0; j <= n - 1; j++)
        {
            q[j] = -a[n - j - 1] / a[n];
        }
        for (j = n; j <= n * n - 1; j++)
        {
            q[j] = 0.0;
        }
        for (i = 0; i <= n - 2; i++)
        {
            q[(i + 1) * n + i] = 1.0;
        }
        i = hhqr(q, n, xr, xi, eps);

        return (i);
    }

    /*
    // 多项式方程求根QR方法例
      int main()
      { 
          int i, n;
          double xr[6],xi[6],eps;
          double a[7]={-30.0,10.5,-10.5,1.5,4.5,-7.5,1.5};
          eps=0.000001;  n=6;
          i=qrrt(a,n,xr,xi,eps);
          if (i>0)
          { for (i=0; i<=5; i++)
              cout <<"x(" <<i <<") = " <<xr[i] <<"   J " <<xi[i] <<endl;
          }
          return 0;
      }
    */
}

